import numpy as np
import pandas as pd

import scipy.optimize
import scipy.stats as st

import iqplot

import tqdm

import bokeh.io
import bokeh.plotting
bokeh.io.output_notebook()

Loading BokehJS ...

df = pd.read_csv('../data/singer_transcript_counts.csv', comment='#')

df.head()

	Rex1	Rest	Nanog	Prdm14
0	11	34	39	0
1	172	91	33	5
2	261	70	68	0
3	178	54	88	1
4	129	54	41	0

n = df['Nanog'].values

p = iqplot.ecdf(
    n,
    x_axis_label='nanog count'
)

bokeh.io.show(p)

def log_like_nbinom(params, n):
    """Log likelihood for repeated measurements from NBinom."""
    alpha, b = params

    if alpha <= 0 or b <=0:
        return -np.inf

    beta = 1 / b
    return np.sum(st.nbinom.logpmf(n, alpha, beta/(1+beta)))


def neg_log_like_nbinom(params, n):
    return -log_like_nbinom(params, n)


def nbinom_mle(n):
    nbar = np.mean(n)
    s2 = np.var(n)
    beta0 = nbar / (s2 - nbar)
    alpha0 = nbar * beta0

    params0 = np.array([alpha0, 1 / beta0])

    res = scipy.optimize.minimize(
        neg_log_like_nbinom, params0, args=(n,), method='Powell'
    )

    alpha, b = res.x

    return alpha, b


def gen_nbimon(alpha, b, size):
    beta = 1 / b
    return np.random.negative_binomial(alpha, beta / (1 + beta), size=size)


def draw_parametric_bs_reps_mle(
    mle_fun, gen_fun, data, args=(), size=1, progress_bar=False
):
    params = mle_fun(data, *args)

    if progress_bar:
        iterator = tqdm.tqdm(range(size))
    else:
        iterator = range(size)

    return np.array(
        [mle_fun(gen_fun(*params, size=len(data)), *args) for _ in iterator]
    )

bs_reps = draw_parametric_bs_reps_mle(
    nbinom_mle,
    gen_nbimon,
    n,
    size=10_000,
    progress_bar=True
)

100%|████████████████████████████████████████████████████████| 10000/10000 [00:47<00:00, 208.62it/s]

np.percentile(bs_reps[:,0], [2.5, 97.5])

array([1.07966995, 1.49956886])

np.percentile(bs_reps[:,1], [2.5, 97.5])

array([56.40912929, 84.05777092])

import bebi103

rg = np.random

def gen_nbimon(params, size, rg):
    alpha, b = params
    beta = 1 / b
    return rg.negative_binomial(alpha, beta / (1 + beta), size=size)

bs_reps = bebi103.bootstrap.draw_bs_reps_mle(
    nbinom_mle,
    gen_nbimon,
    n,
    size=10_000,
)

 12%|███████                                                  | 1246/10000 [00:20<01:01, 142.34it/s]

bs_reps

array([[ 1.35241872, 61.79183402],
       [ 1.37863287, 52.93115143],
       [ 1.15899676, 74.42410383],
       ...,
       [ 1.20938007, 64.19263231],
       [ 1.33370964, 70.62928845],
       [ 1.32687391, 62.26093127]])

df_res = pd.DataFrame(data=bs_reps, columns=['α*', 'β*'])

p = bebi103.viz.corner(
    samples=df_res,
    parameters=['α*', 'β*'],
    plot_ecdf=True
    # show_contours=True,
    # levels=[0.95],
)

bokeh.io.show(p)

p = iqplot.ecdf(
    n,
    x_axis_label='nanog count',
    conf_int=True
)

beta = 1 / b

n_theor = np.arange(0, 501)
cdf_theor = st.nbinom.cdf(n_theor, alpha, beta/(1+beta))
p.circle(n_theor, cdf_theor, color='orange')

bokeh.io.show(p)

np.random.negative_binomial(alpha, beta/(1+beta), size=len(n))

array([157,  78, 151,  99,  82, 108, 130, 136,  36,  78, 186,  64, 137,
        44,  18, 136,  82, 180,  47,  65,  62,  28,   9,  49,  77, 230,
        25,  54, 125,  65,  28,  71,  55, 107,  18, 199,  20,  40, 153,
         0,  16, 140,  28,  85,  84, 223, 195, 293,  33, 172, 138,  47,
        13,  18,  62,  11,  52,  10, 175,  29,   9,  36, 146,  69, 104,
        73, 144,  67, 134,  16,  89,   7,  27,  76,  32,  59,  12,  77,
       126, 197,  56,  12,  44,  37, 190,  78,  49,  73,  54,   7, 133,
        61,  96, 350,  41,  43,   4, 121,  81, 136,  62,  28, 131,  44,
       169,  62,  37, 261,  41,  95,  33,   9, 142,  64, 129,  68,   1,
        32,  95,  96,  57,  23, 199,  33,  57, 125, 104,  39, 179,   5,
        56,  18,  95,  46,  46,  77, 133,  81, 165, 162, 122,  35,  80,
        32,  72,  27,   7,  84, 104,  18,  56, 171,  28, 213,  99,  42,
        65,  97,  53,  16,   9,  17,  33,  76,  94,  57, 133,  53,  46,
       124,  24,  60, 205, 150, 107,  76,  76,  50,   4, 334,  75,  82,
        14, 100, 132,  58,  30,  19,  24,  81,  33,  38, 120,  94, 129,
        38,  73,  24, 117,  15,  66,  23,   7, 190, 140,  59, 164,  41,
       156,   9,  53,   9,  29, 105,  45,  46, 129,   8, 300,   6,  27,
        10,  34, 155,  69,   5,  28,  20, 150, 131,  56,  76, 204,  45,
       171,  60,  20, 121, 107, 213, 218,  46,  53,  13, 138,  29, 161,
       104, 350,  10, 185,  16,  86,  48,  30,  60, 159,  76, 263,  59,
       209,  85, 166,   2, 164,  22,  38,  68, 221,  64,  40,   6, 101,
       227,  73, 126,  24, 161,  27])